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We provide an efficient method for computing the maximum likelihood mixed quantum state (with 
density matrix p) given a set of measurement outcome in a complete orthonormal operator basis 
subject to Gaussian noise. Our method works by first changing basis yielding a candidate density 
matrix p which may have nonphysical (negative) eigenvalues, and then finding the nearest physical 
state under the 2-norm. Our algorithm takes at worst 0(d 4 ) for the basis change plus 0(d 3 ) for 
finding p where d is the dimension of the quantum state. In the special case where the measurement 
basis is strings of Pauli operators, the basis change takes only 0(d 3 ) as well. The workhorse of the 
algorithm is a new linear-time method for finding the closest probability distribution (in Euclidean 
distance) to a set of real numbers summing to one. 
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The race is not always to the swift nor the battle to 
the strong, but that's the way to bet — Damon Runyon 

As scientists, we are faced again and again with the 
problem of determining from imperfect data what "re- 
ally happened" in an experiment. Generally, the more 
data one has the better the reconstruction of the true 
events. Even so, the view remains imperfect, so one typ- 
ically tries only to determine what was the event most 
likely to have led to the observed data. When quantum 
mechanics is considered, the situation is harder on the 
experimentalist, since even the results of perfectly per- 
formed measurements may have probabilistic results. 

Nevertheless, one can still determine the quantum 
state of a system by performing many experiments on 
identically prepared systems and building up good statis- 
tics on the outcomes. If the set of experiments is infor- 
mationally complete then the mixed state density matrix 
describing the system can be determined. This is called 
quantum state tomography [TH3]- A completely deter- 
mination of the state would require an infinite number 
of perfect measurements, which are unobtainable, so in- 
stead we concentrate on finding the maximum likelihood 
state (cf. g]). 

We consider an informationally complete set of mea- 
surements, each performed many times on identically 
prepared quantum systems. From the experimental out- 
comes, we would like to determine the quantum state 
that gives the observed results with highest probability. 
This can be a computationally intensive task. For the 
two qubit experiments of [5], conventional ML solving 
took more time than the experiments themselves. And it 
has been reported in [5] that the ML reconstruction for 8 
qubits in [7] took weeks of computation. Our main result 
is a fast algorithm for reconstructing this state when the 
noise is Gaussian. For 8 qubits our algorithm runs in 
seconds. 

The rest of the paper is organized as follows: First, 
we show that the ML state reconstruction problem with 



Gaussian noise is equivalent to a least-squares minimiza- 
tion problem on quantum states. Next, we prove that the 
minimum takes a particularly simple form. Finally, we 
give a fast algorithm that finds this minimum explicitly 
and benchmark it against direct minimization. 

Reduction to density matrix minimization — Observ- 
ables in quantum mechanics are Hermitian operators 
with thhe expectation value of a Hermitian operator a 
applied to a mixed state p given by Tr(ap). We can rep- 
resent the result of an imperfect measurement of such 
an expectation subject to additive Gaussian noise with 
variance v as a probability density function p{m\p) = 

_J_ p -[m-Tr(*p)] 2 /(2v) 
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Given an orthonormal Hermitian operator basis 
{o~i}f =1 on d x d matrices (with Tr[<7j<7j] = dSij), and 
a particular set of measured values my corresponding to 
jth measurement result of expectation value o~i applied 
to the "true state" po, we want to find the mixed state p, 
a trace 1 Hermitian matrix with only nonnegative eigen- 
values, maximizing the likelihood function 
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Here n is the number of measurement results for the ex- 
pectation of <7j and m, = Y^j=i m ij/ n ^ s the average 
value of those results. The same p that maximizes C will 
minimize the log likelihood function 



^log = }\rnj - Tr ((Tip)] 



(3) 



Working in the operator basis of the {<Ji}s is not con- 
venient, but fortunately the distance is just the Hilbcrt- 
Schmidt, or 2-norm, which is basis independent. We 
show this here for completeness: 
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Lemma: J2i( m i 
and n = Tr[aip\. 
Proof: 



\fi— pW^/d where = Tr[eri/x] 



Im-pII! = Tr (M-p) 2 

= Tr [ i m i ~ Tr P a i )°if 
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= ^Tr[(mi - Tr pai)o-i(mj - Tr pa 

ij 

= )Xrnj - Tr pai)(rrij - Tr pa j) Tr[aiaj] 

ij 

= d (to,; — Tr poi ) 2 

□ 

The matrix p — (1/d) TOjiTj can be thought of as 
the experimentally noisy view of the density matrix po- 
Note that it is trace one by construction, but may have 
negative eigenvalues. Calculation of p from the to^s is 
a change of operator basis, and in general requires time 
0(d 4 ) (there are d 2 values of i and each <7j is a d x d 
matrix). This will actually be the limiting step in our 
overall algorithm, as all other steps with be 0(d 3 ) or 
better. In many cases of interest, however, the operator 
basis change can be done more quickly. This will be 
true whenever the matrices representing the ct^s in the 
canonical basis are sparse. In particular, if the ct^s are 
tensor products of the Pauli matrices 



{a ,a 1: a 2 ,a 3 } = { 



1 
1 



1 

1 
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} 



on n qubits, so that d = 2™, each er^ has only d nonzero 
elements and the change of basis can be carried out in 
0(d 3 ) steps. 

Now, after the change of basis from the m^s to p, 
our original maximum-likelihood problem has been trans- 
formed into the following: 

Subproblem 1: Given a trace-one Hermitian matrix p, 
find the closest density matrix p (a trace-one Hermitian 
matrix with only nonnegative eigenvalues) under the 2- 
norm: 
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This is immediately familiar as a least squares min- 
imization problem, for which standard minimizer pack- 
ages are well suited. Indeed, this is how the problem is 
often solved in practice. Unfortunately, finding the solu- 
tion can be computationally intensive. In standard state 
reconstruction algorithms, this is the most expensive step 
by far (see Fig. [2]) . 

Simple form for the minimum — To improve matters, 
since the 2- norm is basis independent, we can work in 
the eigenbasis of p. Then we observe that the optimum 



p is diagonal in this basis as well. This is immediate from 
the form of ([5| , where any off diagonal terms could only 
contribute positive amounts to the sum. Thus, the prob- 
lem reduces to finding the eigenvectors and eigenvalues 
of p and picking the d non-negative eigenvalues for p that 
minimize Eigensystem solutions are 0(d 3 ) and good 
packages exist [5]. 

We are left with a minimization over d variables, ef- 
fectively the square root of the difficulty of the original 
problem. Call the eigenvalues of p,p Pi,\ and arrange 
them such that pi > pi+i- We now want to minimize 
Si(Aj - Pi) 2 such that J^i Pi = J2i ^ = 1 and ^ > °- 
Using the method of Lagrange multipliers to impose this 
constraint, and writing Xi = x 2 to enforce non-negativity 
of Xi we write the objective function 



Differentiating with respect to Xi we have 



<9A 



A(x 2 - p i )x l - 2Lxi = 



(6) 
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This equation has two solutions, either Xi = or 

x 2 = L/2 + pi . (8) 

Note that L doesn't depend on i so each Xi = x 2 is either 
set to zero or given by pi plus the very same number. To 
evaluate L we must pick a set s = {i such that Xi =/= 0}. 
Then, summing (j8| over i we have 1 = J2i x \ = |s|L/2 + 
E ies Pi or 



£ / 2 = h£. 
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and 
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Lemma: Consider an i and j with p t > pj with i £ s 
and j £ s. Then Aj, the distance function for this case, 
is always less than or equal to Aj, the distance function 
for a new set s' = s + {j} — {i}. 
Proof: We can write 



A; 



and for the case with j g s and i ^ s 
/■'.■ - P., > 2 

If we had A,- < Aj this would imply 



A; = |s| (^^-X-pj+pi + j: 
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pj - pi - {sKpi + pj) 
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Then we would have 

L Mi - Mi 



A, 



< 



(|s| - l)(Mj - Mi) 



2|s| 



< (14) 



because fa > Mr But Aj must be nonnegative, therefore 
A is never decreased by moving i into and j out of s. 
□ 

The lemma tells us that all the A^s that are zero are 
together at the end, matching up with the smallest fas. 
Thus, rather than the 2 d possible choices for s we need 
only decide where to put the break between zero and 
nonzero AjS, for which there are only d choices. 

Next, we show that the choice of s should be the largest 
set satisfying the constraint that all the A^s are nonnega- 
tive. Starting from Eq. (11), we imagine removing some 



element j from from s. Then 



A' = 

and 

A' - A,; 
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In other words, setting any more of the A^s to zero than 
necessary increases the distance function. We are now 
ready to give an algorithm for Subproblem 1. 
Fast algorithm for Subproblem 1 — 

1. Calculate the eigenvalues and eigenvectors of fa 
Arrange the eigenvalues is order from largest to 
smallest. Call these fa, \fa), 1 < i < d. 

2. Let i = d and set an accumulator a = 0. 

3. If Mi + a/i is non-negative, go on to step 4. Oth- 
erwise, set Xi = and add Mi to a. Reduce i by 1 
and repeat step 3. 

4. Set Xj = fij + a/i for all j < i. 

5. Construct p = J2i Ai|Ai)(Aj|. 

Fig. [I] works through an example of this algorithm. 

Efficiency of the algorithm — The slowest step is step 
1, solving the eigensystem, which is 0(d 3 ) (the standard 
libraries such as LAPACK start and are limited by re- 
ducing a Hermitian matrix to tridiagonal form using the 
Householder method which is 0(d 3 ) 9]). Step 2 is obvi- 
ously constant time, and step 3 and 4 together are easily 
seen to be 0(d). Step 5 involves a choice of whether 
one wants the answer in the eigenbasis of fa in which 
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FIG. 1: Example of our algorithm for subproblem 1: (a) We 
start with fa = 3/5, fa = 1/2, /X3 = 7/20, (14 = 1/10, 
— —11/20, and accumulator a — 0. (b) Since ^5 + a/5 
is negative, A5 is set to and a to —11/20. (c) Since /m + a/4 
is negative, A4 is set to and a = —9/20. (d) Finally, Ai, A2, 
and A3 each have a/3 = —3/20 added to the corresponding 
fi. The final result is Ai = 9/20, A 2 = 7/20, A3 = 1/5, and 
A 4 = A 5 = 0. 



case it has already been computed, or in some other ba- 
sis, in which case it is 0(d 3 ) or better [10]. Thus, the 
overall complexity is 0(d 3 ). The actual run-time of an 
implementation will depend primarily on the eigensys- 
tem solver used. Fig. [2] compares the run time of our 
algorithm with that of a traditional ML optimization. 

Discussion — Our argument that the problem of ML 
state reconstruction reduces to finding the closest density 
matrix under the 2-norm to a nonphysical candidate ma- 
trix ju. depended on the measurements of the expectation 
values of Tr (o~ip) in Eq. ([T]) having the same Gaussian 
variance for all i. In practice, this may not be the case. 
One solution to this is to equalize the variances by sim- 
ply performing more of the noisier measurements. For 
example, in circuit quantum electrodynamics (supercon- 
ducting qubits) quantum non-demolition measurement of 
a Paulis are performed with v = 1/ (Tt) where T is the 
measurement rate and r is the measurement time. 

We have focused only on state reconstruction, but it 
is often desired to perform process tomography |llj . that 
is, to determine the quantum input /output relation (a 
trace-preserving completely-positive map) implemented 
by an apparatus. Due to the Choi-Jamiolkowski isomor- 
phism [12l[T3] such a map can be represented by a density 
matrix. Thus, the generalization of our result to process 
tomography should be straightforward. 

Finally, we note a connection to classical probability 
theory. If one considers the eigenvalues fa (some of which 
may be negative) as a noisy view of a probability distri- 
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FIG. 2: Run time for reconstruction of random n-qubit pure 
states mixed with the identity and subjected to Gaussian 
noise on Pauli measurements. Given a set of measurement 
outcomes, we generate a candidate matrix fj, consistent with 
the data but with potentially negative eigenvalues. We use 
two techniques to generate the maximum likelihood density 
matrix p. The top, dashed line, is Matlab's fminsearch used 
to minimize Tr[Qu — p) 2 ] directly. The lowest, dotted line, is 
our algorithm for Subproblem 1. The middle, solid line, is the 
Subproblem 1 algorithm together with the basis change from 
the measurement outcomes to p. All timings were performed 
in Matlab on a single core of an Intel E8400 CPU running at 
3GHz. 

bution, then the algorithm starting from step 2 is an al- 
gorithm for finding the nearest proper probability distri- 
bution. Furthermore, if the noise is Gaussian, this finds 
the maximum likelihood probability distribution. 
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